package analysis;


/**
 * @author Mathieu Rousseau
 * @email mathieu.rousseau3@mail.mcgill.ca
 */


import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Iterator;
import utilities.PDButility;


/** ------------------------------------------------------------------------------------------------------------------
 * start of StructureAnalyzer class
 * ------------------------------------------------------------------------------------------------------------------  **/
public class StructureAnalyzer {
	
	
	/*  global parameters  */
	private String inputPath = "StructureAnalyzer/Input/";
	private String outputPath = "StructureAnalyzer/Output/";
	private double[][] IFs;
	private ArrayList<Integer> featureLocations;
	private ArrayList<String> featureLabels;
	private ArrayList<Integer> genomicMapping;
	private double[][] structure;
	private double[][] pairwiseDistances;
	private ArrayList<double[][]> poolPairwiseDistanceBetweenFeatures;
	private PDButility pdb = new PDButility();
	/*  end of global parameters  */

	
	/*  start of constructor method  */
	public StructureAnalyzer(){
		//  empty constructor
	}/*  end of constructor method  */
	
	
	/*  start of constructor method  */
	public StructureAnalyzer(String featureFilename){
		
		setFeaturesFromFile(featureFilename);
		
	}/*  end of constructor method  */
	
	
	/*  start of setInputPath method  */
	public void setInputPath(String inputPath){
		this.inputPath = inputPath;
	}/*  end of setInputPath method  */
	
	
	/*  start of getInputPath method  */
	public String getInputPath(){
		return inputPath;
	}/*  end of getInputPath method  */

	
	/*  start of setOutputPath method  */
	public void setOutputPath(String outputPath){
		this.outputPath = outputPath;
	}/*  end of setOutputPath method  */
	
	
	/*  start of getOutputPath method  */
	public String getOutputPath(){
		return outputPath;
	}/*  end of getOutputPath method  */
	
	
	/*  start of setIFsFromFile method  */
	public void setIFsFromFile(String IFsFilename){
		//  find the min and max indices
		int minIndex = Integer.MAX_VALUE;
		int maxIndex = Integer.MIN_VALUE;
		
		String line = null;
		try{
			BufferedReader br = new BufferedReader(new FileReader(inputPath + IFsFilename));			
			while((line = br.readLine()) != null){
				String[] values = line.split("\t");
				if(values.length != 5){
					System.out.println("Error: the IFs file does not have exactly five (5) entries per line.\nSystem.exit(0)");
					System.exit(0);
				}
//				System.out.println(line);
				int index1 = Integer.parseInt(values[0]);
				int index2 = Integer.parseInt(values[1]);
				if(index1 < minIndex){
					minIndex = index1;
				}
				if(index2 < minIndex){
					minIndex = index2;
				}
				if(index1 > maxIndex){
					maxIndex = index1;
				}
				if(index2 > maxIndex){
					maxIndex = index2;
				}
			}
			br.close();
		}catch (IOException e) {
			e.printStackTrace();
			System.out.println("Error: cannot find IFs file with filename \"" + IFsFilename + "\"\n" +
			"System.exit(0)");
			System.exit(0);
		}
//		System.out.println("Opened and read file: " + IFsFilename);
		
//		System.out.println("The minimum index is: " + minIndex);
//		System.out.println("The maximum index is: " + maxIndex);
		
		//  make the counting start from zero (0) (same as for the indices of the points in the structure)
		IFs = new double[maxIndex][maxIndex];
		
		line = null;
		try{
			BufferedReader br = new BufferedReader(new FileReader(inputPath + IFsFilename));			
			while((line = br.readLine()) != null){
				String[] values = line.split("\t");
				if(values.length != 5){
					System.out.println("Error: the IFs file does not have exactly five (5) entries per line.\nSystem.exit(0)");
					System.exit(0);
				}
//				System.out.println(line);
				int index1 = Integer.parseInt(values[0]) - 1;		//  make the indices start at zero (0)
				int index2 = Integer.parseInt(values[1]) - 1;		//  make the indices start at zero (0)
				IFs[index1][index2] = Double.parseDouble(values[2]);
//				System.out.println(index1 + "-" + index2 + ":\t" + IFs[index1][index2]);
			}
			br.close();
		}catch (IOException e) {
			e.printStackTrace();
			System.out.println("Error: cannot find IFs file with filename \"" + IFsFilename + "\"\n" +
			"System.exit(0)");
			System.exit(0);
		}
		System.out.println("Opened and read file: " + IFsFilename);		
	}/*  end of setIFsFromFile method  */


	/*  start of setIFs method  */
	public void setIFs(double[][] IFsParam){
		IFs = new double[IFsParam.length][IFsParam[0].length];
		
		for(int i=0; i<IFsParam.length; i++){
			for(int j=0; j<IFsParam[0].length; j++){
				IFs[i][j] = IFsParam[i][j];
			}
		}
	}/*  end of setIFs method  */

	
	/*  start of getIFs method  */
	public double[][] getIFs(){
		return IFs;
	}/*  end of getIFs method  */
	
	
	/*  start of setFeaturesFromFile method  */
	public void setFeaturesFromFile(String featureFilename){
		featureLocations = new ArrayList<Integer>();
		featureLabels = new ArrayList<String>();
		
		String line = null;
		try{
			BufferedReader br = new BufferedReader(new FileReader(inputPath + featureFilename));			
			while((line = br.readLine()) != null){
				String[] values = line.split("\t");
				if(values.length != 2){
					System.out.println("Error: the features files does not have exactly two (2) entries per line.\nSystem.exit(0)");
					System.exit(0);
				}
//				System.out.println(line);
				featureLocations.add(Integer.parseInt(values[1]));
				featureLabels.add(values[0]);
			}
			br.close();
		}catch (IOException e) {
			e.printStackTrace();
			System.out.println("Error: cannot find feature file with filename \"" + featureFilename + "\"\n" +
			"System.exit(0)");
			System.exit(0);
		}
		System.out.println("Opened and read file: " + featureFilename);
		
		featureLocations.trimToSize();
		featureLabels.trimToSize();
	}/*  end of setFeaturesFromFile method  */

	
	/*  start of setFeatureList method  */
	public void setFeatureList(ArrayList<Integer> featureList, ArrayList<String> labelList){
		featureLocations = new ArrayList<Integer>();
		featureLabels = new ArrayList<String>();
		for(Integer location: featureList){
			featureLocations.add(location);
//			System.out.println(location);
		}
		featureLocations.trimToSize();
		for(String label: labelList){
			featureLabels.add(label);
//			System.out.println(label);
		}
		featureLabels.trimToSize();
	}/*  end of setFeatureList method  */
	
	
	/*  start of getFeatureList method  */
	public ArrayList<Integer> getFeatureList(){
		return featureLocations;
	}/*  end of getFeatureList method  */
	
	
	/*  start of getFeatureLabels method  */
	public ArrayList<String> getFeatureLabels(){
		return featureLabels;
	}/*  end of getFeatureLabels method  */
	
	
	/*  start of setGenomicMappingFromFile method  */
	public void setGenomicMappingFromFile(String genomicMappingFilename){
		genomicMapping = new ArrayList<Integer>();
		
		String line = null;
		try{
			BufferedReader br = new BufferedReader(new FileReader(inputPath + genomicMappingFilename));			
			while((line = br.readLine()) != null){
				String[] values = line.split("\t");
				//  there should be two entries, the first is the atom index in the structure, the second is the genomic index
				if(values.length != 2){
					System.out.println("Error: not exaclty two (2) entries per row in genomic mappings file.\nSystem.exit(0)");
					System.exit(0);
				}
//				System.out.println(line);
				genomicMapping.add(Integer.parseInt(values[1]));
			}
			br.close();
		}catch (IOException e) {
			e.printStackTrace();
			System.out.println("Error: cannot find genomic mapping file with filename \"" + genomicMappingFilename + "\"\n" +
			"System.exit(0)");
			System.exit(0);
		}
		System.out.println("Opened and read file: " + genomicMappingFilename);
		
		genomicMapping.trimToSize();
	}/*  end of setGenomicMappingFromFile method  */
	
	
	/*  start of setGenomicMappingList method  */
	public void setGenomicMappingList(ArrayList<Integer> genomicMappingList){
		genomicMapping = new ArrayList<Integer>();
		for(Integer location: genomicMappingList){
			genomicMapping.add(location);
//			System.out.println(location);
		}
		genomicMapping.trimToSize();
	}/*  end of setGenomicMappingList method  */
	
	
	/*  start of getGenomicMappingList method  */
	public ArrayList<Integer> getGenomicMappingList(){
		return genomicMapping;
	}/*  end of getGenomicMappingList method  */
	
	
	/*  start of setStructureFromFile method  */
	public void setStructureFromFile(String structureFilename){
		structure = pdb.getStructure(inputPath + structureFilename);
		System.out.println("Opened and read file: " + structureFilename);
	}/*  end of setStructureFromFile method  */
	
	
	/*  start of setStructure method  */
	public void setStructure(double[][] structure){
		this.structure = new double[structure.length][3];
		for(int i=0; i<structure.length; i++){
			this.structure[i][0] = structure[i][0];
			this.structure[i][1] = structure[i][1];
			this.structure[i][2] = structure[i][2];
		}
	}/*  end of setStructure method  */
	
	
	/*  start of getStructure method  */
	public double[][] getStructure(){
		return structure;
	}/*  end of getStructure method  */
	
	
	/*  start of findCoordinatesOfGenomicPosition method  */
	public double[] findCoordinatesOfGenomicPosition(int genomicPosition){
		//  NOTE: this method assumes that the genomic mappings are in monotonically increasing order
		
		double[] coordinates = new double[3];
		
		//  make sure that the genomic position is within the range of the structure
		if((genomicPosition < Collections.min(genomicMapping)) || (genomicPosition > Collections.max(genomicMapping))){
			System.out.println("Error: the genomic position " + genomicPosition + " is not contained within the structure.");
			System.out.println("System.exit(0)");
			System.exit(0);
		}
		else{
//			System.out.println("The genomic position " + genomicPosition + " is contained within the structure.");
		}
		
		//  find the endpoints of the fragment that contains the genomic position
		if(genomicMapping.size() < 2){
			System.out.println("Error: there is less than one fragment in the structure.\nSystem.exit(0)");
			System.exit(0);
		}
		int fragmentStartIndex = 0;
		Iterator<Integer> genomicMappingIterator = genomicMapping.iterator();
		int startPosition = genomicMappingIterator.next();
		int endPosition = genomicMappingIterator.next();
		while(genomicMappingIterator.hasNext()){
			if((genomicPosition >= startPosition) && (genomicPosition <= endPosition)){
				break;
			}
			else{
				startPosition = endPosition;
				endPosition = genomicMappingIterator.next();
				fragmentStartIndex++;
			}
		}
		//  at this point the genomic position MUST be contained within the region defined by [startPosition,endPosition]
		double sizeOfFragment = endPosition - startPosition;
		double distanceIntoFragment = genomicPosition - startPosition;
		double fractionIntoFragment = distanceIntoFragment / sizeOfFragment;
		
		int fragmentEndIndex = fragmentStartIndex + 1;
		
		double xVector = structure[fragmentEndIndex][0] - structure[fragmentStartIndex][0];
		double yVector = structure[fragmentEndIndex][1] - structure[fragmentStartIndex][1];
		double zVector = structure[fragmentEndIndex][2] - structure[fragmentStartIndex][2];
		
		coordinates[0] = structure[fragmentStartIndex][0] + fractionIntoFragment * xVector;
		coordinates[1] = structure[fragmentStartIndex][1] + fractionIntoFragment * yVector;
		coordinates[2] = structure[fragmentStartIndex][2] + fractionIntoFragment * zVector;
		
		return coordinates;
	}/*  end of findCoordinatesOfGenomicPosition method  */
	

	/*  start of getEuclideanDistanceBetweenGenomicPositions method  */
	public double getEuclideanDistanceBetweenGenomicPositions(int positionOne, int positionTwo){
		double[] P1 = findCoordinatesOfGenomicPosition(positionOne);
		double[] P2 = findCoordinatesOfGenomicPosition(positionTwo);
		
		return getDistanceBetweenCartesianPositions(P1, P2);
	}/*  end of getEuclideanDistanceBetweenGenomicPositions method  */
	
	
	/*  start of getDistanceBetweenCartesianPositions method  */
	private double getDistanceBetweenCartesianPositions(double[] P1, double[] P2){
		double distanceBetweenPoints = 0;
		
		double term = P2[0] - P1[0];		//  x2 - x1
		term *= term;						// (x2 - x1)^2
		distanceBetweenPoints += term;
		term = P2[1] - P1[1];				//  y2 - y1
		term *= term;						// (y2 - y1)^2
		distanceBetweenPoints += term;
		term = P2[2] - P1[2];				//  z2 - z1
		term *= term;						// (z2 - z1)^2
		distanceBetweenPoints += term;
		
		return distanceBetweenPoints;
	}/*  end of getDistanceBetweenCartesianPositions method  */
	
	
	/*  start of getFeatureCoordinatesInStructure method  */
	public double[][] getFeatureCoordinatesInStructure(){
		double[][] featureCoordinates = new double[featureLocations.size()][3];
		for(int i=0; i<featureLocations.size(); i++){
			featureCoordinates[i] = findCoordinatesOfGenomicPosition(featureLocations.get(i));
		}
		return featureCoordinates;
	}/*  end of getFeatureCoordinatesInStructure method  */
	
	
	/*  start of outputFeatureCoordinatesInStructure method  */
	public void outputFeatureCoordinatesInStructure(String filename){
		double[][] featureCoordinates = getFeatureCoordinatesInStructure();
		
		int locationOfExtension = filename.lastIndexOf('.');
		String pdbFilename = filename.substring(0, locationOfExtension) + ".pdb";
		
		pdb.writeAnnotation(outputPath + pdbFilename, featureCoordinates, featureLabels, "annotation");
		
		try{						
			BufferedWriter bw = new BufferedWriter(new FileWriter(outputPath + filename));
			//  write the header line
			bw.write("feature\tx\ty\tz");
			bw.newLine();
			bw.flush();
			for(int i=0; i<featureLocations.size(); i++){
				bw.write(featureLocations.get(i) + "\t" + featureCoordinates[i][0] + "\t" + featureCoordinates[i][1] + "\t" + featureCoordinates[i][2]);
				bw.newLine();
				bw.flush();
			}
			bw.flush();
			bw.close();
		}catch (IOException e) {
			System.out.println("Error: could not output file: " + filename);
			e.printStackTrace();
		}//end of try-catch block
		System.out.println("Output file: " + outputPath + filename);
	}/*  outputFeatureCoordinatesInStructure  */
	
	
	/*  start of calculatePairwiseDistanceBetweenFeatures method  */
	public double[][] calculatePairwiseDistanceBetweenFeatures(){
		//  calculate all of the pairwise distances
		pairwiseDistances = new double[featureLocations.size()][featureLocations.size()];
		for(int i=0; i<featureLocations.size()-1; i++){
			for(int j=i+1; j<featureLocations.size(); j++){
				pairwiseDistances[i][j] = getEuclideanDistanceBetweenGenomicPositions(featureLocations.get(i), featureLocations.get(j));
			}
		}
		return pairwiseDistances;
	}/*  end of calculatePairwiseDistanceBetweenFeatures method  */
	
	
	/*  start of outputPairwiseDistanceBetweenFeaturesAsList method  */
	public void outputPairwiseDistanceBetweenFeaturesAsList(String filename){
		try{						
			BufferedWriter bw = new BufferedWriter(new FileWriter(outputPath + filename));
			//  as a list
			for(int i=0; i<featureLocations.size()-1; i++){
				for(int j=i+1; j<featureLocations.size(); j++){
					bw.write(featureLocations.get(i) + "-" + featureLocations.get(j) + "\t" + pairwiseDistances[i][j]);
					bw.newLine();
					bw.flush();
				}
			}
			bw.flush();
			bw.close();
		}catch (IOException e) {
			System.out.println("Error: could not output file: " + filename);
			e.printStackTrace();
		}//end of try-catch block
		System.out.println("Output file: " + outputPath + filename);
	}/*  end of outputPairwiseDistanceBetweenFeaturesAsList method  */
	

	/*  start of outputPairwiseDistanceBetweenFeaturesAsHeatmap method  */
	public void outputPairwiseDistanceBetweenFeaturesAsHeatmap(String filename){
		try{						
			BufferedWriter bw = new BufferedWriter(new FileWriter(outputPath + filename));
			//  as a heatmap
			//  write the header line
			bw.write("\t");
			for(int i=0; i<featureLocations.size()-1; i++){
				bw.write(featureLocations.get(i) + "\t");
			}
			bw.write(String.valueOf(featureLocations.get(featureLocations.size()-1)));
			bw.newLine();
			bw.flush();
			for(int i=0; i<featureLocations.size(); i++){
				bw.write(featureLocations.get(i) + "\t");
				for(int j=0; j<i; j++){
					bw.write("\t");
				}
				bw.write("0\t");
				for(int j=i+1; j<featureLocations.size(); j++){
					bw.write(pairwiseDistances[i][j] + "\t");					
				}
				bw.newLine();
				bw.flush();
			}
			bw.flush();
			bw.close();
		}catch (IOException e) {
			System.out.println("Error: could not output file: " + filename);
			e.printStackTrace();
		}//end of try-catch block
		System.out.println("Output file: " + outputPath + filename);		
	}/*  end of outputPairwiseDistanceBetweenFeaturesAsList method  */
	
	
	/*  start of calculatePoolAnalysisPairwiseDistanceBetweenFeatures method  */
	public void calculatePoolAnalysisPairwiseDistanceBetweenFeatures(int numberOfStructures){
		
		poolPairwiseDistanceBetweenFeatures = new ArrayList<double[][]>(numberOfStructures);
				
		for(int structureIndex=1; structureIndex<=numberOfStructures; structureIndex++){
			//  get the structure filename (default, structure numbering starting from index one (1))
			String filenameStructure = pdb.getDefaultFilenameForIndex(structureIndex);
			
			//  set the current structure of the analyzer object
			setStructureFromFile(filenameStructure);
			
			//  calculate the pairwise distances between the features for the current structure and add it to the list
			poolPairwiseDistanceBetweenFeatures.add(calculatePairwiseDistanceBetweenFeatures());
		}		
		poolPairwiseDistanceBetweenFeatures.trimToSize();
	}/*  end of calculatePoolAnalysisPairwiseDistanceBetweenFeatures method  */
	
	
	/*  start of outputPoolAnalysisPairwiseDistanceBetweenFeatures method  */
	public void outputPoolAnalysisPairwiseDistanceBetweenFeatures(String filename){
		//  calculate the average and the standard deviation for the pairwise distances across all of the structures in the pool
		double[][] averagePairwiseDistance = new double[featureLocations.size()][featureLocations.size()];
		for(double[][] currentArray: poolPairwiseDistanceBetweenFeatures){
			for(int i=0; i<featureLocations.size()-1; i++){
				for(int j=i+1; j<featureLocations.size(); j++){
					averagePairwiseDistance[i][j] += currentArray[i][j];
				}
			}
		}
		//  if there is more than one (1) structure then calculate the average pairwise distance in the pool
		if(poolPairwiseDistanceBetweenFeatures.size() > 1){
			for(int i=0; i<featureLocations.size()-1; i++){
				for(int j=i+1; j<featureLocations.size(); j++){
					averagePairwiseDistance[i][j] /= (double)poolPairwiseDistanceBetweenFeatures.size();
				}
			}
		}
		else if(poolPairwiseDistanceBetweenFeatures.size() == 1){			//  there is only one structure so it is the average
			averagePairwiseDistance = poolPairwiseDistanceBetweenFeatures.get(0);
		}
		else{
			System.out.println("Error: there are no structures to calculate the average from.\nSystem.exit(0)");
			System.exit(0);
		}
//		System.out.println("-----     AVERAGES     -----");
//		for(int i=0; i<featureLocations.size()-1; i++){
//			for(int j=1; j<featureLocations.size(); j++){
//				System.out.print(averagePairwiseDistance[i][j] + "\t");
//			}
//			System.out.println();
//		}
		
		
		//  only calculate the standard deviation if there are at least three (3) structures in the pool
		double[][] stdevPairwiseDistance = new double[featureLocations.size()][featureLocations.size()];
		if(poolPairwiseDistanceBetweenFeatures.size() >= 3){
			for(double[][] currentArray: poolPairwiseDistanceBetweenFeatures){
				for(int i=0; i<featureLocations.size()-1; i++){
					for(int j=i+1; j<featureLocations.size(); j++){
						double term = currentArray[i][j] - averagePairwiseDistance[i][j];
						term *= term;
						stdevPairwiseDistance[i][j] += term;
					}
				}
			}
			for(int i=0; i<featureLocations.size()-1; i++){
				for(int j=i+1; j<featureLocations.size(); j++){
					stdevPairwiseDistance[i][j] /= (double)(poolPairwiseDistanceBetweenFeatures.size() - 1.0);
					stdevPairwiseDistance[i][j] = Math.sqrt(stdevPairwiseDistance[i][j]);
				}
			}
		}
//		System.out.println("-----     STDEVS     -----");
//		for(int i=0; i<featureLocations.size()-1; i++){
//			for(int j=1; j<featureLocations.size(); j++){
//				System.out.print(stdevPairwiseDistance[i][j] + "\t");
//			}
//			System.out.println();
//		}
		
		
		//  output the average pairwise distances to file
		String outputFilenameAverages = "pairwiseDistancesPool_average_" + filename;
		try{
			BufferedWriter bw = new BufferedWriter(new FileWriter(outputPath + outputFilenameAverages));
			//  as a list
			for(int i=0; i<featureLocations.size()-1; i++){
				for(int j=i+1; j<featureLocations.size(); j++){
					bw.write(featureLocations.get(i) + "-" + featureLocations.get(j) + "\t" + averagePairwiseDistance[i][j]);
					bw.newLine();
					bw.flush();
				}
			}
			bw.flush();
			bw.close();
		}catch (IOException e) {
			System.out.println("Error: could not output file: " + outputFilenameAverages);
			e.printStackTrace();
		}//end of try-catch block
		System.out.println("Output file: " + outputPath + outputFilenameAverages);
		
		
		//  output the stdev pairwise distances to file
		String outputFilenameStdevs = "pairwiseDistancesPool_stdev_" + filename;
		try{
			BufferedWriter bw = new BufferedWriter(new FileWriter(outputPath + outputFilenameStdevs));
			//  as a list
			for(int i=0; i<featureLocations.size()-1; i++){
				for(int j=i+1; j<featureLocations.size(); j++){
					bw.write(featureLocations.get(i) + "-" + featureLocations.get(j) + "\t" + stdevPairwiseDistance[i][j]);
					bw.newLine();
					bw.flush();
				}
			}
			bw.flush();
			bw.close();
		}catch (IOException e) {
			System.out.println("Error: could not output file: " + outputFilenameStdevs);
			e.printStackTrace();
		}//end of try-catch block
		System.out.println("Output file: " + outputPath + outputFilenameStdevs);
	}/*  end of outputPoolAnalysisPairwiseDistanceBetweenFeatures method  */
	
	
	/*  start of compactionProfile method  */
	public void compactionProfile(String filename, int diagonalSizeToExclude){
		try{						
			BufferedWriter bw = new BufferedWriter(new FileWriter(outputPath + filename));
			//  write the header line
			bw.write("index1\tindex2\tEuclid\tIF\tExcluded diagonal of size=" + diagonalSizeToExclude);
			bw.newLine();
			bw.flush();
			//  go through every pair of points in the structure and output the matched IF and Euclidean distance
			for(int i=0; i<structure.length-1; i++){
				for(int j=i+1; j<structure.length; j++){
					
					if((j-i) > diagonalSizeToExclude){
						double IF = IFs[i][j];
						double Euclid = getDistanceBetweenCartesianPositions(structure[i], structure[j]);
					
						bw.write((i+1) + "\t" + (j+1) + "\t" + Euclid + "\t" + IF);		//  output the indices with counting starting at one (1) and not zero (0)
						bw.newLine();
						bw.flush();
					
//						System.out.println(i + "-" + j + ":");
//						System.out.println("\tEuclid=" + Euclid);
//						System.out.println("\tIF=" + IF);					
					}
				}
			}
			bw.flush();
			bw.close();
		}catch (IOException e) {
			System.out.println("Error: could not output file: " + filename);
			e.printStackTrace();
		}//end of try-catch block
		System.out.println("Output file: " + outputPath + filename);		
	}/*  end of compactionProfile method  */
	
	
}/*  end of StructureAnalyzer class  */
